using System;
using System.Collections.Generic;
using System.Text;

using Athena.Core;

namespace Athena.Tools
{
    /// <summary>
    /// This class contains a set of static methods that generate a transcription rate
    /// expression when provided with the correct Kd values. The class also provides 
    /// ready-made rate expressions such as AND, OR, NOT, XOR, NAND. These ready-made
    /// rate expressions need a list of nodes and stoichiometry values, and they return
    /// an Edge that contains the neccessary modifiers and parameters. 
    /// </summary>
    public class TranscriptionRate
    {
        public static double[] RPolBindingProb; //2^n
        public static double[,] ProteinBindingProb; //n x 2^n
        public static int[,] booleanTable;
        public static string[] ProteinNames;
        public static List<double> ParameterValues = new List<double>();

        /// <summary>
        /// Give the species default names - s1,s2,s3...
        /// These names will be used in the rate expression
        /// </summary>
        /// <param name="N"></param>
        public static void DefaultNames(int N)
        {
            ProteinNames = new string[N];
            for (int i = 0; i < N; ++i)
                ProteinNames[i] = string.Format("s{0}", i+1);
        }

        /// <summary>
        /// initialize the rate-building arrays such as species names and number
        /// of configurations, etc.
        /// </summary>
        /// <param name="nodes">a array of nodes</param>
        /// <param name="stoichiometry">an array of corresponding stoichiometry values</param>
        public static void InitializeWithNodes(Node[] nodes, double[] stoichiometry)
        {
            if (nodes == null) return;

            List<string> names = new List<string>();
            for (int i=0; i < nodes.Length; ++i)
            {
                if (stoichiometry != null && stoichiometry.Length == nodes.Length && stoichiometry[i] != 1.0d)
                    names.Add(nodes[i].ID + "^" + stoichiometry[i]);
                else
                    names.Add(nodes[i].ID);                
            }

            TranscriptionRate.FillBooleanTable(names.Count);
            for (int i = 0; i < names.Count && i < ProteinNames.Length; ++i)
            {
                ProteinNames[i] = names[i];
            }
        }
        /// <summary>
        /// initialize the rate-building arrays such as species names and number
        /// of configurations, etc.
        /// </summary>
        /// <param name="edge">an Edge</param>
        public static void InitializeWithEdge(Edge edge)
        {
            if (edge == null) return;

            List<string> names = new List<string>();

            foreach (Node n in edge.Reactants)
            {
                double stoic = edge.GetReactantStoichiometry(n);
                if (stoic != 0)
                {
                    if (stoic != 1.0d)
                        names.Add(n.ID + "^" + stoic);
                    else
                        names.Add(n.ID);
                }
            }
            foreach (Node n in edge.Modifiers)
            {
                double stoic = edge.GetModifierStoichiometry(n);
                if (stoic != 0)
                {
                    if (stoic != 1.0d)
                        names.Add(n.ID + "^" + stoic);
                    else
                        names.Add(n.ID);
                }
            }
            
            TranscriptionRate.FillBooleanTable(names.Count);
            for (int i = 0; i < names.Count && i < ProteinNames.Length; ++i)
            {
                ProteinNames[i] = names[i];
            }
        }

        /// <summary>
        /// initialize the array of possible configurations
        /// </summary>
        /// <param name="N">the number of species</param>
        public static void FillBooleanTable(int N)
        {
            booleanTable = new int[(int)Math.Pow(2, N), N];  //a table will all possible bound/not-bound conformations
            FillBooleanTable(new int[N], 0, 0);   //make that table
            if (ProteinNames == null || ProteinNames.Length != N)    //default names
            {
                DefaultNames(N);
            }
            ParameterValues.Clear();
        }

        /// <summary>
        ///fills booleanTable with 1/0 in increasing (binary numbers) order
        /// </summary>
        /// <param name="prot"></param>
        /// <param name="depth"></param>
        /// <param name="k"></param>
        /// <returns></returns>
        public static int FillBooleanTable(int[] prot, int depth, int k)
        {
            if (depth >= prot.Length)
            {
                for (int i = 0; i < prot.Length; ++i)
                    booleanTable[k, i] = prot[i];
                return k + 1;
            }
            else
            {
                prot[depth] = 0;
                k = FillBooleanTable(prot, depth + 1, k);
                prot[depth] = 1;
                k = FillBooleanTable(prot, depth + 1, k);
            }
            return k;
        }

        /// <summary>
        /// converts a binary array into a number
        /// </summary>
        /// <param name="nums"></param>
        /// <returns></returns>
        public static int ArrayToInt(int[] nums)
        {
            int x = 0;
            for (int i = 0; i < nums.Length; ++i)
                if (nums[i] == 1)
                    x += (int)Math.Pow(2, i);
            return x;
        }

        /// <summary>
        /// determines the prob of a set of proteins being bound to the node
        /// </summary>
        /// <param name="prot"></param>
        /// <returns></returns>
        public static double ProteinBindingCoeff(int[] prot)  //prot = set of protein (binary array)
        {
            int k = 0, sum = 0;
            for (int i = 0; i < prot.Length; ++i) //check if only one protein is bound
            {
                if (prot[i] == 1) k = i;
                sum += prot[i];
            }
            if (sum == 1) return ProteinBindingProb[k, 0];  //then prob = [p]*[G]*k

            double answer = 0;
            for (int i = 0; i < prot.Length; ++i)   //do every combination of p_i + G(others) -> G(all p)
                if (prot[i] == 1)
                {
                    int[] temp = (int[])prot.Clone();
                    temp[i] = 0;
                    int x = ArrayToInt(temp);

                    answer += ProteinBindingProb[i, x] * ProteinBindingCoeff(temp);

                    break;
                }
            return answer;
        }

        /// <summary>
        /// checks that all paths leading to this configuration have same rates
        /// </summary>
        /// <param name="prot"></param>
        /// <returns></returns>
        public static double AllPathsEquivalent(int[] prot)
        {
            int k = 0, sum = 0;
            for (int i = 0; i < prot.Length; ++i) //check if only one protein is bound
            {
                if (prot[i] == 1) k = i;
                sum += prot[i];
            }
            if (sum == 1)
                return ProteinBindingProb[k, 0];  //base case

            double answer = 0, s = 0, s0 = 0;
            for (int i = 0; i < prot.Length; ++i)   //do every combination of p_i + G(others) -> G(all p)
                if (prot[i] == 1)
                {
                    int[] temp = (int[])prot.Clone();
                    temp[i] = 0;
                    int x = ArrayToInt(temp);
                    if (ProteinBindingProb[i, x] > 0)
                    {

                        s0 = s;
                        s = AllPathsEquivalent(temp) * ProteinBindingProb[i, x];
                        if (s < 0 || (s0 > 0 && s != s0))
                        {
                            return -1;
                        }
                        answer += s;
                    }
                }
            return answer;
        }

        /// <summary>
        /// Obtain the string version of ProteinBinding
        /// </summary>
        /// <param name="prot"></param>
        /// <param name="bindingProb"></param>
        /// <returns></returns>
        public static string sProteinBinding(int[] prot, string bindingProb)
        {
            string answer = bindingProb;
            if (answer == "0") return "0";
            for (int i = 0; i < prot.Length; ++i)   //do every combination of p_i + G(others) -> G(all p)
                if (prot[i] == 1)
                {
                    answer += string.Format("*{0}", ProteinNames[i]);
                }
            return answer;
        }

        /// <summary>
        /// get the equation for the activity of a gene
        /// </summary>
        /// <param name="N">number of species involved</param>
        /// <returns></returns>
        public static string sGetProb(int N) //n is the number of promoters
        {
            if (RpolEntirelyZero()) return "0";

            FillBooleanTable(N);

            StringBuilder denominator = new StringBuilder("("), numerator = new StringBuilder("(");

            int n = booleanTable.GetLength(1), numRows = booleanTable.GetLength(0);   //total possible configurations of bound promoters

            for (int i = 0; i < numRows; ++i) //for each R + G(p1p2...) -> RpolG(p1p2...)
            {
                StringBuilder s1 = new StringBuilder(), s2 = new StringBuilder(); //temp

                int[] temp = new int[n];
                for (int j = 0; j < n; ++j) temp[n - 1 - j] = booleanTable[i, j]; //copy each row

                string bindingProb = sProteinBinding(temp, "k" + (ParameterValues.Count + 1));  //get the prob. for that row as string   

                if (AllPathsEquivalent(temp) < 0)
                    return "Model disobeys thermodynamics!";

                if (bindingProb.Length > 1 || i == 0)
                {
                    if (bindingProb.Length > 2)
                    {
                        ParameterValues.Add(ProteinBindingCoeff(temp)); //get the prob. for that row
                        s1.Append(bindingProb); //a*p1
                        //if (RPolBindingProb[i] > 0) s1.Append("*");
                    }
                    /*
                    if (RPolBindingProb[i] > 0)                        
                        if (RPolBindingProb[i] == 1)
                            s1.Append("(1 + RPol)");
                        else
                        {
                            ParameterValues.Add(RPolBindingProb[i]);
                            s1.Append(string.Format("(1 + {0}*RPol)", "k" + ParameterValues.Count));  //a*p1*(1 + b*Rpol)
                        }
                         
                    else*/
                    if (i == 0)
                        {
                            s1.Append(string.Format("1 "));
                        }
                }


                if (RPolBindingProb[i] > 0) //i.e. trascriptionally active
                {
                    if (i == 0)
                    {
                        //ParameterValues.Add(RPolBindingProb[i]);
                        //s2.Append("k" + ParameterValues.Count);
                        s2.Append(" 1 ");
                    }
                    else
                    if (bindingProb.Length > 1)
                    {
                        s2.Append(bindingProb);
                        //s2.Append("*");
                    }
                    
                        /*
                        if (bindingProb.Length > 1 || i == 0)
                        if (bindingProb.Equals("k" + ParameterValues.Count))
                            s2.Append("RPol");
                        else
                            s2.Append(string.Format("{0}*RPol", "k" + ParameterValues.Count));  //active conformations in the numerator
                         */
                }
                if (s1.Length > 1)
                {
                    if (denominator.Length > 1) denominator.Append(" + ");
                    denominator.Append(s1);
                }
                if (s2.Length > 1)
                {
                    if (numerator.Length > 1) numerator.Append(" + ");
                    numerator.Append(s2);
                }

            }
            numerator.Append(")/ ");
            denominator.Append(")");
            numerator.Append(denominator);

            string str = numerator.ToString();

            ParameterValues.Add(10d);
            str = ("k" + ParameterValues.Count) + "*" + str;

            return str;
        }

        /// <summary>
        /// presence of Rna polymerase
        /// </summary>
        /// <returns></returns>
        private static bool RpolEntirelyZero()
        {
            for (int i = 0; i < RPolBindingProb.Length; ++i)
            {
                if (RPolBindingProb[i] > 0)
                    return false;
            }
            return true;
        }

        /// <summary>
        /// Returns the rate if all the TFs are additive (upregulate)
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string LinearUpRegulation(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                RPolBindingProb[i] = temp;
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns the rate if all the TFs are additive (downregulate)
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string LinearDownRegulation(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                RPolBindingProb[i] = 1/(1 + temp);
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns the rate for OR control
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string OR(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                if (temp > 0)
                    RPolBindingProb[i] = 10;
                else
                    RPolBindingProb[i] = 0;
            }

            return sGetProb(nBindingSites);
        }

        /// <summary>
        /// Returns the rate for AND control
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string AND(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                if (temp >= nBindingSites)
                    RPolBindingProb[i] = 10;
                else
                    RPolBindingProb[i] = 0;
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns the rate for NOR control
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string NOR(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                if (temp > 0)
                    RPolBindingProb[i] = 0;
                else
                    RPolBindingProb[i] = 10;
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns the rate for NAND control
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string NAND(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                if (temp < nBindingSites)
                    RPolBindingProb[i] = 10;
                else
                    RPolBindingProb[i] = 0;
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns the rate for XOR control
        /// </summary>
        /// <param name="nBindingSites"></param>
        /// <returns></returns>
        public static string XOR(int nBindingSites)
        {
            RPolBindingProb = new double[(int)Math.Pow(2, nBindingSites)];
            ProteinBindingProb = new double[nBindingSites, (int)Math.Pow(2, nBindingSites)];
            FillBooleanTable(nBindingSites);

            for (int i = 0; i < booleanTable.GetLength(0); ++i)
            {
                int temp = 0;
                for (int j = 0; j < nBindingSites; ++j)   //get the vector for the row
                {
                    temp += booleanTable[i, j];
                    ProteinBindingProb[j, i] = 1; //set the value
                }

                if (temp == 1)
                    RPolBindingProb[i] = 10;
                else
                    RPolBindingProb[i] = 0;
            }

            return sGetProb(nBindingSites);
        }
        /// <summary>
        /// Returns an edge that represents the OR operation on the given nodes. stoichiometry can be null
        /// </summary>
        /// <param name="nodes"></param>
        /// <param name="stoichiometry"></param>
        /// <returns></returns>
        public static string OR(Node[] nodes, double[] stoichiometry)
        {
            InitializeWithNodes(nodes, stoichiometry);
            string rate = OR(nodes.Length);
            return rate;
        }
        /// <summary>
        /// Returns an edge that represents the AND operation on the given nodes. stoichiometry can be null
        /// </summary>
        /// <param name="nodes"></param>
        /// <param name="stoichiometry"></param>
        /// <returns></returns>
        public static string AND(Node[] nodes, double[] stoichiometry)
        {
            InitializeWithNodes(nodes, stoichiometry);
            string rate = AND(nodes.Length);
            return rate;
        }
        /// <summary>
        /// Returns an edge that represents the NAND operation on the given nodes. stoichiometry can be null
        /// </summary>
        /// <param name="nodes"></param>
        /// <param name="stoichiometry"></param>
        /// <returns></returns>
        public static string NAND(Node[] nodes, double[] stoichiometry)
        {
            InitializeWithNodes(nodes, stoichiometry);
            string rate = NAND(nodes.Length);
            return rate;
        }
        /// <summary>
        /// Returns an edge that represents the NOR operation on the given nodes. stoichiometry can be null
        /// </summary>
        /// <param name="nodes"></param>
        /// <param name="stoichiometry"></param>
        /// <returns></returns>
        public static string NOR(Node[] nodes, double[] stoichiometry)
        {
            InitializeWithNodes(nodes, stoichiometry);
            string rate = NOR(nodes.Length);
            return rate;
        }
        /// <summary>
        /// Returns an edge that represents the XOR operation on the given nodes. stoichiometry can be null
        /// </summary>
        /// <param name="nodes"></param>
        /// <param name="stoichiometry"></param>
        /// <returns></returns>
        public static string XOR(Node[] nodes, double[] stoichiometry)
        {
            InitializeWithNodes(nodes, stoichiometry);
            string rate = XOR(nodes.Length);
            return rate;
        }

        /// <summary>
        /// Copies the parameter values generated by the transcription rate tool to the
        /// edge. The edge must have a ParentModule in order for this method to do anything
        /// </summary>
        /// <param name="edge"></param>
        public static void CopyParameterValues(Edge edge)
        {
            if (edge == null) return;

            Module module = edge.ParentModule;
            if (module == null) return;

            List<string> parms = edge.ParameterNames;
            parms.Sort();

            if (module != null)
                for (int i = 0; i < parms.Count && i < Tools.TranscriptionRate.ParameterValues.Count; ++i)
                {
                    module.SetParameter(parms[i], Tools.TranscriptionRate.ParameterValues[i]);
                }
        }
    }
}

